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Abstract 

We study the fluctuational dynamics of a tagged base-pair in double stranded DNA. We cal- 
culate the drift force which acts on the tagged base-pair using a potential model that describes 
interactions at base pairs level and use it to construct a Fokker-Planck equation. The calculated 
displacement autocorrelation function is found to be in very good agreement with the experimental 
result of Altan-Bonnet et. al. Phys. Rev. Lett. 90, 138101 (2003) over the entire time range of 
measurement. We calculate the most probable displacements which predominately contribute to 
the autocorrelation function and the half-time history of these displacements. 

PACS numbers: 87.15.-v, 05.40.-a, 87.10.-he,02.50.-r 
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DNA double stranded helical structure is stabilized by the hydrozen bonding between 
complementary bases and the stacking between neighbouring bases [1]. In physiological 
solvent conditions the average value of these interactions for a base pair is of the order 
of few ksT (thermal energy) [2] and thermal fluctuations can lead to local and transitory 
unzipping of the double strands [3,4] . The co-operative opening of a sequence of consecutive 
base pairs leads to formation of local denaturation zones (bubbles). As an AT base pair 
connected by two hydrogen bonds needs less energy to unzip compared to a GC base pair 
which is connected by three hydrogen bonds, initiation of a bubble generally takes place 
in an AT rich region. A DNA bubble consists of flexible single stranded DNA and its size 
fluctuates by zipping and unzipping of base pairs at the two zipper forks where the bubble 
connects to the intact double strands. The average size of a bubble depends on the sequence 
of base pairs, temperature and ionic strength and varies from few broken base pairs at room 
temperature to few hundred open base pairs close to melting temperature [5,6]. 

The formations of bubble at room or physiological temperatures are rare and intermittent 
with life times of the order of millisecond [4]. The occurrence of such bubble domains is 
important as the opening of dsDNA base pairs by breaking the hydrogen bonds between 
complementary bases disrupts the helical stack and may initiate biological processes of tran- 
scription, replication and protein binding [7,8]. Prom physics point of view, DNA bubbles 
offer a quasi one -dimensional system for the study of fluctuational dynamics. 

In an experiment by Altan-Bonnet et. al. [4] the dynamics of a single bubble in three 
synthetic DNA constructs having the same GC rich region but different AT base pairs regions 
have been investigated by fluorescence correlation spectroscopy (FCS). In the middle of the 
AT region a T base pair was tagged with a fluorophore while the neighbouring T base of 
the other strand was tagged with a quencher. The correlation spectrum of fluctuating base 
pairs was monitored. The dynamics was found to follow a multi-state relaxation kinetics in 
a wide temperature range with a characteristic time scale in the range of 20 — lOOyUS. Several 
theoretical models [4,9,12,13]have recently been proposed to explain the observed multi-state 
breathing dynamics. In one of these models [9] the bubble free energy that corresponds to 
a bubble of inflnitely large size [10,11] and which accuracy for a bubble of few broken base 
pairs, to best of our knowledge, is not established has been used. Other theoretical models 
include discrete master equations and stochastic Gillespie schemes [4,12,13]. 

In this Letter, we develop a general theory to study the fluctuational dynamics of a tagged 
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base pair by means of a Fokker-Planck equation based on a potential field which acts on the 
base pair and which we obtain by integrating out the degrees of freedom of all base pairs of 
a dsDNA except those associated with the tagged one. We use the simple potential model 
of Peyard-Bishop-Dauxious (PBD) [14] to represent the interactions in dsDNA at base pairs 
level. 

The PBD model reduces the degrees of freedom of DNA to a one -dimensional chain of 
effective atom compounds describing the relative base pair separation t/j from the ground 
state position yi=0. The potential of the model is written as 

i 

+k/2{1 + pe-"^y^+y'-^\yi - (1) 

where N is the number of base pairs, summation on the r.h.s. is over all base pairs of the 
molecule and y^={yi}, the set of relative base pair separations. The first term of Eq.(l) is 
the Morse potential that represents the hydrogen bonds between the bases of the opposite 
strands and the second term represents the stacking interaction between adjacent base pairs. 
The values of parameters found by Campa and Giansanti [15] are k — 
and a = 0.35A~^ for the stacking part, while for the Morse potential Dat = 0.05 eV, qat 
= 4.2A^^ for an AT base pair and Dec = 0.075 eV and acc = 6.9A^^ for a GC base pair. 

We now consider one of the DNA molecules (named A18) investigated by Altan-Bonnet 
et. al. [4] and take the VJth base pair counted from the 5'— end as the tagged base pair. The 
interactions in the molecule is represented by the PBD model. We add a harmonic potential 
Uh(yN)—0(yN — '^)cy% where c ~ 1.0 x 10~^A~^ and 6(y) a Heaviside step function at the 
terminal GG base pair to avoid the complete separation of the two strands. In experiment 
[4] this was achieved by attaching a hairpin loop of 4T. The potential felt by the tagged 
base pair at a separation y from the ground state y = is found from the relation 

V{y) = -kBT[lnZr,{yo) - /nZ,(0)] (2) 

where = / Ug^dyiSiy,, - y)exp [-/3C/(y^)] 

Zn{0) ^inl,dyiS{yn - 0)exp [-/3C/(?/^)] 
are the constrained partition function integrals, S is the Dirac function and /3 = (kBT)^^ ■ 
For the PBD model the calculation of a partition function integral reduces to multipli- 
cation of N matrices. The discretization of the co-ordinate variable and introduction of a 
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proper cut off on the maximum values of y's determines the size of the matrices. We have 
taken — 2A and 120A as the lower and upper limit of integration for each co-ordinate variable 
and discretized space using the Gaussian- Legendre method with number of grid points equal 
to 900. Note that the values of the partition function integrals and therefore the values of 
V{y) are independent of the hmit of integration. We show in Fig.l the value of V{y) as a 
function of y at 45°C. At a separation y the base pair feels a drift force F — —dV{y)/dy to- 
wards the origin y = 0. As shown in the inset of Fig.l this force has a minimum a.t y = 0.2 A. 
This minimum corresponds to a force barrier which has been observed in theoretical investi- 
gation of force induced unzipping of a dsDNA in the constant extension ensemble [16,17,18] 
and is attributed to a combination of the force needed to break the hydrogen bonds and the 
force needed to overcome the entropic barrier of the stacking interaction [16]. 
The dynamics of the base pair may be described by the Langevin equation 

ft = + ^ ' < > (^) = '^^ksT5{t) (3) 

where F is a transport coefficient of dimension time/mass and F/cfiTa^y of dimension 
1/time. Eq.(3) describes a one-dimensional random walk in a potential V{y). We use a at 
and rkBTa\rp make, respectively, distance and time dimensionless. The Fokker-Planck 
equation corresponding to (3) is found to be 



dP _ d 
dt dy 



-dpV{y) 
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9 P 



dy 

where P{y, yo] t) is the probability density of the random walkers. 

We assume that if separation y reduces to zero at time t', it will not contribute to 
autocorrelation function defined as C{t) =< y(t)y{0) > — < y >^ for t > t' and similarly 
any new fluctuational opening which appear after t = will not contribute to C{t). Thus 
for purposes of computing the autocorrelation function we place an absorbing wall at y = 
0, i.e.P{y — 0,t) — 0. In addition to this we may require P{y — L,t) — 0, where L depends 
on the size of the dsDNA molecule or on any other condition which limits the size of the 
bubble. The problem of calculating the autocorrelation function C{t) therefore reduces to 
finding how many walkers of an ensemble of random walkers distributed according to thermal 
equilibrium distribution at t = are still present at time t and have not been absorbed by 
the wall at y = [19]. 
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When a substitution P = exp[—j3V{y)/2\ip is used Eq.(4) reduces to 

^2 
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This is the imaginary time Schrodinger equation for a particle of mass 1/2 in the potential 
v{y). Let (f)m{y) denote the eigenfunctions of the operator H, Hcf)^^ — with ^^(y = 

0) = and / dy4>^{y)(f)m'{y) — 5mm' Then expanding ip{y,t) in terms of eigenfunctions (pm 
and using the initial condition P{y,yo,t = 0) = S{y — yo) the transition probability from 
initial separation y^ to a final separation y at time t is found to be 

y{y) + v{yoy 



P{y:yo:t) = exp 



-(3- 



m 



-Emt 



0m(l/)0m(l/o) 



(7) 



For initial distribution of separation yQ we choose the Boltzmann factor B{y) — 
Aexp{—(3V{y)) where A = 1/ dyoexp{—/3V{yQ)) is a normalization factor. If we start 
with the equilibrium distribution function at time t = 0, the distribution function at time 
t,P{y,t) is 

y{y) + V{yoy 



P{y,t) = A^^e-^-' f dyo 



exp 



-(5- 



(l>m{y)(l)*m{yo) 



(8) 



P{y, t) measures the survival probability. 
For the autocorrelation function we get 



C{t)= f"^ P{y,t)dy = Aj^e-""-' C 
Jo tr^ ^0 



-PViy)/2^ 



i{y)dy 



(9) 



The values of (t)m{y) and Em of the operator H in Eq.(5) are determined numerically 
using a method developed by Sethia et.al.[20]. As shown in Fig. 2(a), v{y) is attractive at 
small y, rises to a (repulsive) maximum at y = 0.2A and then decays to zero as y increases. 
The maximum in v{y) corresponds to the minimum in F shown in Fig.l. For small values of 
m , (f)miy) remains confined (see Fig. 2(b)) in a region of separation which values are smaller 
than the length of the molecule L. After the first three eigenvalues which values are Eq — 
0.0028, El = 0.0080, E2 = 0.0125, Em is found to increase with AEm = Em+i -Em^ 0.0043 
for m < 50. The free particle in a box like behaviour is found only after m > 100 and 
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therefore the values of 0m and E.^ depend on the value of L only after m > 100. We have 
varied L from 80A — 120A and found that the values of C{t) and P{y, t) do not change. The 
values given in Fig.3 and Fig.4 correspond to L = lOOA which approximately measures the 
length of the dsDNA molecule of 29 base pair. 

In Fig.3 the rescaled autocorrelation function g{u) — G{t/ti/2) where G{t) — C{t)/C{0) 
and ti/2 is such that G(ti/2) = 0.5 [4], is plotted as a function of rescaled time t/ti/2. When 
this figure is compared with the one given in [4] we find a very good agreement over the 
entire time range of measurement. If we choose TkBTa^jp = 10^s~^ and plot G{t) as a 
function of time t in ms the resulting curve is found to be in very good agreement with the 
corresponding curve given in [4]. 

In Fig.4(a) we show the distribution function P(y, t) which gives the probability of sep- 
aration y of the tagged base pair at time t. From the figure we find that the most probable 
separation is y* ^ lA, although the term "most probable" makes less and less sense because 
the peak gets broader and broader. Thus, initially as well as presently small separation 
of the order of lA make the most contributions to the autocorrelation function C{t) at all 
times. This can be understood from the nature of the drift force F{y) (shown in Fig.l) which 
favours small separation. Since small separations have larger Boltzmann weights initially, 
they dominate C{t) at all times. If we plot Iny* vs Int we find a straight line having a slope 
equal to 1/6. Thus the most probable displacements of the base pair depends on time as 

The half time history of a random walker that is at at t = and at t is defined as [19] 



H{y, t/2\y\ t; y\ t = 0) = P{y\ t\y, t/2)P{y, t/2; y\ t = 0) 
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We plot the half time distribution as a function of y for t = 5, 10 and 20 in Fig 4(b). While 
values of y* corresponding to these times are 0.84, 1.01 and 1.24 the peak in H are found 
respectively at 0.90, 1.12 and 1.48 which are somewhat larger than the corresponding values 
of y*.The half width of the distribution H{y,t/2\y*,t;y*,0) is found to be narrower than 
that of P{y,t). Therefore the most probable way for a displacements of size y* formed at 
t = to survive until a time t is that they first grow larger than y* and then shrink back to 
the original size. 

In conclusion; we developed a theory to describe the multi-state relaxation dynamics of 
a tagged base pair of dsDNA. We used a potential model which describes interactions in 



dsDNA at base pairs level and calculated the drift force which acts on the base pair and 
drives it to its equilibrium position. The dynamics is governed by the Langevin equation 
with Gaussian white noise. We derived the associated Fokker-Planck equation and with 
suitable transformation reduced it to an imaginary time Schrodinger equation for a particle 
of mass 1/2. We found the eigenvalues and eigenfunctions of the operator using a numerical 
method described in [20] . The calculated displacement autocorrelation function is found to 
agree with experimental result for the entire time range of measurement. The most probable 
displacements which contribute predominately to short as well as long times are found to 
be small, of the order of lA. The half time distribution of these displacements which show 
how the most probable displacements behave between time t — and t are calculated. 
The method developed here is equally applicable to homogeneous and heterogeneous DNA 
molecules. 
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FIG. 1: The effective potential felt by the tagged base pair in a dsDNA molecule (A18 of [4]) of 
29 base pairs at separation y at 45°C In the inset the drift force F{y) = —dV{y)/dy which drives 
the base pair to the equilibrium position is plotted. 
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FIG. 2: (a)Tlic potential v{y) of Eq.(6) at 45°C. The (repulsive) maximum in v{y) correspond 
to the minimum in the drift force F(y). (b) Results of first few eigenfunctions as a function of y 
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FIG. 3: Rescaled autocorrelation function g{u) = G{t/ti/2) where G{t) = C(i)/C(0) and ti/2 is 
such that G{ti/2) = 0.5 as a function of i/ii/2 at 45^(7. These notations are same as used in [4]. 
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FIG. 4: (a)Results for the distribution P{y, t) as a function of separation y at 45°C for several time 
t which are expressed in unit of {TkBTa\j,)~^ . The peak in P{y,t) represents the most probable 
separations, (b) Results for the half-time distribution H(y,t/2\y* ,t;y* ,t = 0) as a function of 
separation y for i = 5, 10, 20 for which y* = 0.84, 1.01, 1.24. The peak in H is found at separation 
larger than the corresponding value of y*. 
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